Eighth-order phase-field-crystal model for two-dimensional crystallization 



O 

> 
O 

CO 
CN 



A. Jaatinen 1 ' 2 and T. Ala-Nissila 1 ' 3 
1 Department of Applied Physics, Aalto University School of Science, P.O. Box 11000, FI-00076 Aalto, Finland 

2 Department of Materials Science and Engineering, 
Aalto University School of Science, P.O. Box 16200, FI-00076 Aalto, Finland and 
3 Department of Physics, Brown University, Providence RI 02912-1843 
(Dated: 19 November 2010) 

We present a derivation of the recently proposed eighth order phase field crystal model [Jaatinen 
et al. , Phys. Rev. E 80, 031602 (2009)] for the crystallization of a solid from an undercooled 
melt. The model is used to study the planar growth of a two dimensional hexagonal crystal, and the 
results are compared against similar results from dynamical density functional theory of Marconi 
and Tarazona, as well as other phase field crystal models. We find that among the phase field 
crystal models studied, the eighth order fitting scheme gives results in good agreement with the 
density functional theory for both static and dynamic properties, suggesting it is an accurate and 
computationally efficient approximation to the density functional theory. 
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I. INTRODUCTION 

Understanding crystal formation from an undercooled 
melt is of significant academic interest due to the com- 
plex phenomena involved in crystallization, and also of 
practical interest due to its relevance to a vast amount 
of industrial processes. During the past decade, rapidly 
evolving progress in microscopic understanding of phe- 
nomena involved in solidification has followed the intro- 
duction of the phase field crystal (PFC) model [1,2]. This 
model was first introduced as a phenomenological exten- 
sion of the traditional phase field models Q such that 
the order parameter field exhibits the crystalline nature 
of the underlying crystal lattice. The most significant 
advantage of PFC type of models over the traditional 
phase field models is that including the periodic struc- 
ture of the field in the model will naturally result in in- 
clusion of many crystal structure related properties, such 
as elasticity, plasticity and grain boundaries [l| . Since its 
introduction, the PFC model has been applied to mod- 
eling elastic and plastic deformation of materials 0, 0] , 
dislocation dynamics @ , crystal growth @, 0] , static and 
dynamic properties of driven two dimensional overlayers 

Because of the periodic nature of the order parameter 
field in the PFC model, it is not hard to come up with an 
intuitive interpretation that the field must be related to 
the atomic number density of the underlying system. On 
the other hand, studies of classical density functional the- 
ory (DFT), most commonly in the context of inhomoge- 
nous liquids [ll[ , have aimed at a microscopic derivation 
of the static (and more recently dynamic [Hj]) proper- 
ties of the systems under study by using the microscopic 
density as a field variable in the theory. The extension 
of this approach to crystallization is known as the DFT 
of freezing, which has, in its many forms, been applied 
to study freezing of many different classical systems with 
a varying degree of success [H, [l4|. In 2007, Elder et 
al. [15j introduced the idea that assuming the field under 



study in the PFC model to be linearly proportional to the 
atomic number density in the DFT, the PFC model can 
be viewed as a simplified version of the DFT, and showed 
that the free energy functional used in PFC studies can 
be derived from the DFT by making certain approxi- 
mations. Wu and Karma [161 ] introduced another way 
of obtaining the parameters for the PFC model using a 
DFT-like approach. In a recent paper, the stren gths and 
weaknesses of the approaches proposed in Refs. [15| and 
[l6| were studied, and a new variant of the PFC model 
known as the eighth-order fitting model (EOF), which 
reproduces certain thermodynamic properties of the ma- 
terial under study significantly more accurately than the 
previously proposed methods, was proposed [17| . More 
recently, the EOF model has been applied to study grain 
boundaries and homogenous nucleation [l9[ of body- 
centered cubic iron. 

In the present work, we will present an alternative in- 
terpretation of the EOF model, in which the field under 
study in the EOF is related to the physical atomic num- 
ber density through a convolution, which filters out the 
sub-atomic wavelength Fourier modes from the atomic 
number density. Using this interpretation, we are able to 
derive the free energy of the EOF model in a way which 
we believe is more consistent with the original DFT than 
the previously presented derivations. Predictions of the 
EOF model are then tested against the DFT and other 
related PFC models. While most studies of the PFC 
model's connection to DFT have concentrated on the free 
energy, i.e. the static properties of the model, we will also 
compare the predictions for crystal growth rates, which is 
a dynamical phenomenon. Similar comparison between 
crystal growth rates in DFT and PFC models has previ- 
ously been published by van Teeffelen et al. who found 
the growth rates of colloidal crystals in the early stages of 
solidification agree relatively well between the DFT and 
their choice of PFC models (not including the EOF) Q ■ 
In the present work, we aim at a more thorough assess- 
ment of the crystal front propagation in the DDFT and 
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PFC models. Instead of the initial growth rate, we aim 
at assessing the steady state front propagation velocity in 
both the diffusion controlled and interface kinetics con- 
trolled regimes using both DDFT and EOF models. In 
addition, the results from these models are compared to 
results of a more traditional fourth-order PFC approach, 
and the "PFC1" model that van Teeffelen et al. proposed, 
and argued to be a more accurate approximation to DFT 
than the other model utilized in their comparison [7] . We 
find that among the PFC models studied, the EOF gives 
results in best agreement with the DFT for both static 
properties and crystal growth. 



II. THEORY 

We study solidification dynamics in a two dimensional 
ensemble of Brownian particles interacting via an inverse 
twelfth-power pair potential, 



v(r) = e(~) 



(1) 



where r is the interparticle separation, e sets the energy 
scale and a is the diameter of the particles. Due to scaling 
properties, all structural and thermodynamic properties 
of this model system only depend on the scaled density 



P = 



k B T 



pa 



(2) 



where fcgT is the thermal energy scale and p — N/A 
is the number of particles per unit area. According to 
molecular dynamics simulations of Broughton et al., the 
equilibrium state of the system is a fluid at densities up 
to pi = 0.987, and a hexagonal solid at densities above 
p s = 1.006, while between these two densities the equilib- 
rium state is a coexistence of the solid and liquid phases 
[20T l2l| . As we assume that hydrodynamic interactions 
and the inertial terms can be neglected, the equations of 
motion for the particles are given by 



r 4 = 7 - 1 (F t +f l ), i=l...N, 



(3) 



where the dot denotes time derivative, 7 is a friction co- 
efficient, Fi is the force from the other particles and an 
external field acting on particle i, and is a Gaussian 
random force that fulfills the fluctuation-dissipation the- 



A. Dynamical density functional theory 

In the DDFT approach, instead of solving the positions 
of individual particles as a function of t, one derives an 
equation of motion for the one-particle density defined 
by 



p(r,t)=£>(i 



>(*))) 



(4) 



where the angular brackets denote a noise average 
Marconi and Tarazona [12J have shown that from the 
equations of motion ([3]), one can derive an equation of 
motion for p(r, t) through a coordinate transformation 
and a subsequent noise averaging. In another procedure, 
Archer and Evans [22j have derived the same equation 
of motion by using the Smoluchowski equation as their 
starting point. The equation of motion for p(r, t) result- 
ing from both of these derivations reads 



7 -1 V 



p(r,t)V 



SF[p(r,t)] 
Sp(r,t) 



(5) 



where F [p(r, t)] is the Helmholtz free energy of the sys- 
tem described by a density field p(r,t). As noted in the 
recent work of Ramos et al. this equation of motion can 
also be obtained in the overdamped limit of a more gen- 
eral equation of motion for the number and momentum 
densities, if the effective Hamiltonian is replaced by the 
free energy and thermal fluctuations are ignored L'-jj . 
The free energy F consists of three parts, F — Fid + 

Fex F^ 

contribution 



, ., i -l xs : where the first term represents the ideal gas 



Fid = k B T J dvp{v) (ln(p(r)A|) - l) , 



(6) 



where At is thermal de Broglie wavelength and the sec- 
ond term is a contribution from an external field, 



drp(r)u ex (r), 



(7) 



where u ex (r) is an external field acting on the particles. 
The third part of F is the excess, which is due to the 
interparticle interactions. For this quantity, exact ex- 
pressions only exist for a very limited range of cases, and 
more generally, approximations will have to be made [Tl| . 
In the present work, we will use the simplest possible non- 
local approximation that is an expansion of F xs in powers 
of Ap = p — po around a uniform reference density po, 
where 



F xs =-k B T J drc«(p )Ap(r) 



k R T 



J J dvdv' Ap(r) C ( 2 >(|r-r'|, Po )Ap(r'), 



(8) 



where e- n ' are called n th order direct correlation func- 
tions. The quantity c^ 2 ' is the Ornstein-Zernike direct 
correlation function that can be obtained from exper- 
iments, computer simulations or a number of approxi- 
mate closure relations to the Ornstein-Zernike equation 
1 1| . In the present work, we will utilize the well-known 
Percus-Yevick closure relation with the pair potential Eq. 
(P) to obtain . The reference density po is chosen such 
that F has two equal minima: The trivial uniform mini- 
mum p(r) = po and another where p(r) has a hexagonal 
structure (the external field u ex is set to zero). This is 
the procedure taken in most DFT studies of freezing, and 
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the resulting po is interpreted as the freezing point of the 
liquid. Indeed, the free energy functional defined by Eqs. 
([6]) and fg| is the simplest free energy functional used in 
static DFT studies of freezing, and its success has varied 
from case to case ([14j and references therein). The free 
energy defined by Eqs. (|6]) and dHJ is also the free energy 
functional to which Elder et al. attempted to connect 
the free energy used in PFC studies [15j . 

Putting our free energy functional (with u ex = 0) to- 
gether with Eq. (JSJ) , and rescaling the density field vari- 
able as p(r, t) = po(l + n(r, t)), we end up with an equa- 
tion of motion 



dn ~ 
7T = V 2 n- 

OT 



(l + n)V / dr'Cdr-r'lHr') , (9) 



-(2) 



where the rescaled time r = r y^ 1 ksTt and C = poc y 
Very few studies of the dynamics of crystallization using 
Eq. ([9]) are found in the literature. Van Teeffelen et al. 
have studied the dynamics of colloidal crystal nucleation 
and found results that seemed to agree well with the re- 
sults of molecular dynamics simulations [24| . The same 
group has studied the initial growth velocity of a colloidal 
crystallization front using Eq. © and compared the re- 
sults against the results obtained from the PFC model 
To our knowledge, no other attempts to assess solid- 
ification dynamics by direct application of Eq. @ exist 
in the literature. 



B. Eighth order phase field crystal model 

The eighth order phase field crystal model (EOF) was 
recently presented in Ref. [l7} in the context of quan- 
titative modelling of body-centered cubic iron. In that 
case, the model was shown to reproduce the anisotropic 
solid-liquid interfacial free energies, bulk moduli of solid 
and liquid phases, and the equilibrium coexistence gap 
between them to a quantitatively satisfactory precision. 
While some of these properties had been reproduced in 
previous versions of phase field crystal models, the com- 
bination seemed inaccessible without the eighth order ex- 
tension. [13]. In a subsequent study, it was shown that 
the EOF is also capable of describing grain boundary 
energies of bcc iron quantitatively [l8[ . 

Despite its quantitative success, there remain open is- 
sues in the EOF model. First, the previously presented 
derivation of the EOF is based on the assumption that 
the field n, related to the atomic number density as in 
the case of DFT, is small, such that the logarithmic term 
in the free energy could be expanded as a power series. 
However, it is well known that the actual atomic number 
density in the solid resembles a set of highly localized 
Gaussians, which does not agree with the assumption of 
n being small. Neither does it suggest that the non-local 
part of the free energy could be assumed local in fc-space, 
as must be assumed, when expanding the direct correla- 
tion function in fc-space around its main maximum. The 
equilibrium n-field resulting from these approximations 



is highly localized in fc-space, having little resemblence to 
the Gaussian-like density obtained from the DFT. Most 
prominently, in the supposedly empty spaces between the 
lattice sites, the field n resulting from the EOF reaches 
values smaller than — 1 , corresponding to unphysical neg- 
ative densities. 

In what follows, we will present an alternative deriva- 
tion of the EOF model which, although far from exact, 
avoids the previously mentioned caveats. The key to our 
derivation is the obvious fix in the interpretation of the 
field n used as the field variable in the EOF model: In- 
stead of insisting that the field n in the EOF model would 
be locally and linearly related to p in the same way as 
in the DDFT model, we will assume they are related 
through a weighing function w as 



n(r) 



Po 1 



dr'w(\r-r'\)(p(r)- p ). (10) 



In the current approach, the weighing procedure intro- 
duced in Eq. (fTU| will act as a Fourier filter, cutting off 
the short wavelength modes of p that have very small 
amplitude in the periodic solutions of the PFC model. 
In other words, the field n in the bulk of the solid phase 
will closely resemble the one-mode approximation, 



ra(r) ps n + 



2u cos 



Any 
v3d„, 



27ra;\ / 2ny 
— 2 cos I — — I cos I —= 

dnn J V V 3d nn 



(11) 



where uq is the fractional density change and d nn is the 
nearest-neighbor distance, even though the underlying 
density field were highly peaked around lattice sites. 

A convenient choice for w is a function, whose Fourier 
transform is given by 



m = 



' 1 - C(k) 
l-C EOF {k) 



(12) 



where CEOF(k) is the "approximation" to C(fc) intro- 
duced in [1711. 



CW(*0 = C(k m )-E s 



1-2 



fc 5 



h 2 

n. .. . 



ft.-. 



(13) 

where k m is the position of the main peak in C(fc), Es 
is chosen such that second derivatives at the peak of the 
original and the approximated curves are equal, and Eb 
is then chosen such that the infinite-wavelength (fc = 0) 
limits are equal. The function CEOF{k) will follow the 
original C(k) very closely from the k — limit up to 
the main peak at fc m , after which the two curves diverge, 
C(fc) approaching zero in an oscillatory fashion, while 
CeofOs) falls in the negative infinity, such that w(k) 
defined by Eq. (TT2"|) will fall close to zero rapidly after the 
main peak, providing the desired Fourier filter property 
mentioned earlier. 
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What makes the choice of w(k) defined by Eq. (fT3|) 
particularly convenient is that the linear part of the free 
energy (linearization of Eqs. (O and (jHJ) can now be 
exactly written as 



Po 2 
1 



drdr' Ap(r) (8(r - r') - C(r - r')) Ap(r') 



drn(r) (1 - C BOF (V 2 )) n(r), 



(14) 



where 5 is the Dirac delta function and CeofC^ 2 ) is the 
inverse Fourier transform of CEOF(k), i.e. Eq. (1131) with 
— fc 2 replaced by V 2 . 

Unfortunately for the non-linear parts of the free en- 
ergy, 



(iF nl =/3{F - F u 



= / dr U(r)ln(p(r)/po)-Ap(r)- 



Ap(r) 
2po 



(15) 



the situation is not nearly as trivial, because the expres- 
sions will become both non-local and non-linear. How- 
ever, as we know that Fu n defined by Eq. (TUT) already 
provides us with both a preferred wavelength of fluctua- 
tions in the system, and a large free energy penalty for 
Fourier modes with k >> k m , we argue that if the am- 
plitudes u in Eq. (|11[) vary on length scales larger than 
the range of the weighing function, it may be sufficient 
to aproximate F n i with a functional that is local in terms 
of the field n. To fit this purpose, we postulate a sim- 
ple, local, non-linear functional, consisting of third and 
fourth order terms, 



PFnl,EOF 



Pi) 



dr 



6 



n(r) 3 



12 



n(r) 4 



(16) 



where a and b are phenomenological constants. As a 
consequence of this crude approximation, ignoring prac- 
tically all information about the Fourier modes with 
k >> k m , it is admittedly obvious that the density field 
obtained from a solution of the EOF model through in- 
verting Eq. (1101) will not be an accurate approximation 
to the real underlying density field, unless a more accu- 
rate approximation to F n i will be presented in subsequent 
studies. 

It should be noted that the free energy functional for 
the EOF obtained by summing up Eqs. (jT4"]) and (fTHl) is 
exactly the same as used in previous EOF studies [17l[l8| . 
One crucial difference, in addition to the different approx- 
imations involved, is that the current derivation does not 
suggest that a = b = 1, like the previous version, based 
on Taylor expansion, did. 17] Instead, in the current 
approach, it is clear that the proper way to choose pa- 
rameters a and b is such that free energies of relevant 
density profiles are reproduced as accurately as possible 
(in the end, to correct for the flaws in the derivation, fit- 
ting the parameters a and b with the desired amplitude 



of the solid phase was the approach taken in the previ- 
ous studies as well). In order to gain insight into how 
these parameters should be chosen, consider a density 
field consisting of an infinite set of normalized Gaussians 
in a triangular lattice, 



p( r ) = y - e -«l r - R 



(17) 



where Rj's are positions of lattice points that belong 
in the underlying hexagonal lattice. It is well known 
that in the bulk of the solid, this is a fairly accurate 
approximation to the density profile that results from 
DFT of freezing [l4j]. More specifically, consider the 
case when the density of lattice points 2/(v / 3d 2 m ) = p a 
and the length of the principal reciprocal lattice vector 
|G m | = 47r/(v3rfnTi) coincides with position of the main 
peak in C, i.e. \G m \ = k m . Then, through Eqs. (fTUj) 
and (I12p , the field n will be closely approximated by Eq. 
(fTTI) . with no = 0, and the amplitude u is related to the 
Gaussian parameter a as 



(18) 



It is then straightforward to calculate from Eqs. (fT4|) 
and (fT6|) that free energy of the system described by this 
single-mode n field will be given by 



PFeof 
TV 



3(1 - C(k m ))u 2 - 2au 3 + —bu 



(19) 



In what follows, this expression will be related to the 
small and large a limits of the original free energy, given 
by Eqs. ^ and |8"]). in the Gaussian approximation. 

In the limit of small a, deviations from uniformity in 
the density field are small, and therefore, instead of the 
full free energy functional, it is sufficient to consider the 
linearised version defined by Eq. (TT4]) . 

EL f f drdr f A / x «/ r _ f a _ q, t _ A / A 

Po 2 JJ 



i£(l-C(| G< 



(20) 

where the sum is over all non-zero reciprocal lattice vec- 
tors. Using Eq. (fl8)) . this can be re- written in terms of 
u as 



pF 
Po 



iWl-C(|G,|) 



[Ojj ' 



(21) 



Now it is easily observed that in the limit where u is 
small, the leading contribution to the sum in Eq. (|21j) 
comes from the shortest reciprocal lattice vectors, for 
which | Gi| = |G m |. As the number of vectors in this 
first star of reciprocal lattice vectors is six, we obtain 
exactly the same small-it behaviour from Eqs. (|19p and 
TT\i . This result may not be surprising, due to the exact 
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FIG. 1: Ideal free energy as a fucntion of u in the one-mode 
approximation to the EOF model, as indicated in Eq. (|24|) 
with an = 3/4 and bid = 1/5 (solid line). Dashed lines rep- 
resent the small- and large-ti limits obtained from the origi- 
nal Fid in the Gaussian approximation. For comparison, the 
dashed-dotted line shows the curve expected from a Taylor 
expansion of Fid (i.e. = bid = 1) which is seen to signifi- 
cantly overestimate Fid- 



relation between the linear parts of free energies in the 
models. 

In the limit of large a, i.e. highly localised density 
peaks, the free energy can be accurately evaluated by 
ignoring the overlap between peaks in the ideal gas term 
[14| . resulting in 



(22) 



Using Eq. (fT5)l . this can also be expressed as a function 
of u as 



^ « In -ln(- ln(u))- \ £ (7(|G,|) U 2 ( W. 

(23) 

We note that it might be possible to come up with 
a form of F^eof that would resemble this form more 
closely than the simple fourth order polynomial in Eq. 
(I19p . For the time being, however, we find it sufficient to 
study the implications of Eq. (|23|) to the parameters a 
and b of Eq. (fl6|) . Perhaps the most interesting of these 
implications is the contribution of the ideal gas term in 
the parameters a and b. Ignoring all terms that are pro- 
portional to C(k) in Eqs. ([T9j) and (|23|) leaves us with 
the problem of finding coefficients aid and bid (subscript 
id refers to the ideal free energy) , such that 

3u 2 - 2a ld u 3 + ^-b ld u A » In (-^j - ln(- \n(u)). (24) 

There are, of course, an infinite number of ways to per- 
form this fit. After experimenting with least squares fits 



on different intervals of u, we found that for example by 
choosing aid = 3/4 and bid = 1/5, we obtain a reason- 
ably good agreement between the two curves over a large 
interval of u, as shown in Fig. [TJ Even though details 
in the fitting procedure will affect the obtained numbers 
to some extent, a common feature observed in all reason- 
able fits is that the parameter aid is always on the order 
of 1, while bid is almost an order of magnitude smaller. 
Thus, this argument explains why the parameters that 
were observed to be practicable in our previous studies 
fl7j were so different from unity, suggested by the naive 
Taylor expansion of the log-term. 

Further refinement to the parameters a and b could be 
obtained by studying the terms related to C(k) in Eq. 
(1231) . However, as the expressions for p and n used in 
deriving Eqs. (|19[) and (1231) . we suggest using numerical 
fitting methods in order to achieve maximum accuracy. 
In order to find a functional, that gives such a field n that 
best approximates the original DFT, the most obvious 
choice is to find such a and b that a numerical free energy 
minimization results in the same u as obtained from the 
DFT, and the solid that exhibits this u coexists with the 
liquid phase (given that is the case in the original DFT 
functional as well). 

For dynamics of the EOF model, we use the form 
widely used in PFC studies, i.e. 



— - v 2 ( 5Feof 
dr \ Sn 



1 - C{k n 



n n h — n 

2 3 



Es 



n + E B 



(25) 



where r is defined as earlier. Motivation to choosing 
this equation of motion is that it is probably the mini- 
mum complexity model satisfying the usual requirements 
for conserving the total mass and evolving towards mini- 
mum of the free energy, that also catches approximately 
the same dynamics as the DDFT in the near-uniformity 
limit, for the relevant Fourier modes up to k m (linearized 
version of Eq. © is d T h(k) = — fc 2 (l — C{k))h(k), while 
that of Eq. (|25l) would be the same, but with C replaced 
by Ceof)- For studying solidification, we believe the 
limit of near-uniformity is the dominant factor affecting 
the solidification front velocity, even though some details 
of the dynamics on the solid side of the front may have 
a secondary effect on the front propagation. 



C. Other phase field crystal models 

In addition to comparing the results of the EOF model 
with those of the DDFT model, we will also compare 
their results to two other PFC models presented in the 
literature. The first model we will include is essentially 
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the Swift-Hohenberg (SH) model used in almost all the 
PFC studies up to date. Where the EOF model contains 
gradients up to eighth order in the linear part of the free 
energy, the SH formulation only contains gradients up to 
fourth order. The procedure we use for obtaining param- 
eters for the SH model such that it could be used to model 
real parameters is essentially the one introduced by Wu 
and Karma (l6| . In the framework of the present work, 
we may also view the SH-based approach an approxima- 
tion to the EOF, where Eb — 0, and the parameters a 
and b are fitted through the same procedure as in the 
EOF. The equation of motion for this model, which we 
will call the fourth-order fit (FOF) for the remainder of 
this paper, becomes 

dn 



l-C{km))n--n' 



(26) 



We note here that even though many formulations of the 
SH-based PFC models do not include the third order 
term —an 3 / 6 in the free energy, it has been shown that 
the model without that term can be exactly recovered 
from Eq. (j2"6")l after appropriate scaling of the field vari- 
able and the parameters. [25( 

In addition to the FOF and EOF models, another in- 
teresting PFC model was proposed in the recent work 
of van Teeffelen et al. Q- In the model they call PFC1, 
they start with the DDFT, and approximate the function 
C(k) by expanding it around k m in a fourth order power 
series, in a similar manner as in the FOF model. How- 
ever, as the excess part of the free energy in this approach 
would not be sufficient to stabilize the solid at any rea- 
sonable density 0, Qj} , the excess part of the free energy 
is multiplied by a scaling factor a. Thus, the equation of 
motion in this model becomes 

dn 2 
— = V n- 

OT 



aV 



(1 + n)V <^ C(k m ) - E s 



(27) 

Using arguments based on a single mode approximation 
to the free energy of the solid, van Teeffelen et al. come up 
with a = 1.15 for the case of colloids interacting via r~ 3 
potential which they studied [7|. In the present work, 
we will utilize numerical fitting methods to find an a, 
such that the correct freezing point from the DFT is re- 
produced in the PFC1 model. Van Teeffelen et al. also 
showed that this PFC1 model performs slightly better 
than their "PFC2" model, which is based on FOF with 
a = b = 1 and a = 1.15, in reproducing the initial crys- 
tallization velocities of the DDFT, arguing that the bet- 
ter success is due to fewer approximations made in the 
derivation. Q In the following, we shall see how the 
PFC1 model compares with EOF and FOF for the case 
under study in the present paper. 




FIG. 2: Direct correlation function at pi — 0.9156 as obtained 
from the PY closure relation (solid line), and the fitted EOF 
(dashed line) and FOF (dashed-dotted line) expansions. 



III. STATIC RESULTS AND MODEL 
PARAMETERS 



In order to find the freezing density of the fluid in the 
DFT model, we have calculated the direct correlation 
function for the pair potential defined by Eq. (JlJ at dif- 
ferent densities by using the well-known Ornstein-Zernike 
equation together with the closure relation by Percus and 
Yevick (P Y) [26| . At each density, we then find the non- 
trivial minimum of the free energy, in which the density 
field has a hexagonal structure (at densities where it ex- 
ists), by a similar free energy minimization method as 
utilized in Of the minimization method, it is wor- 

thy of noting that unlike in most DFT studies 0, E3| (and 
like in [13), we do not restrict the calculations to a per- 
fect lattice (i.e. zero vacancy concentration) constraint, 
for the reason that in a dynamical simulation, it is not 
possible to control this issue without modifying the free 
energy. 

Repeating the procedure of finding the direct correla- 
tion function and the solid minimum of free energy at 
many different reference densities, we find that freezing 
occurs at pi = 0.9156, because for that reference density, 
the minimum free energy of the solid equals that of the 
liquid at p = pq. That is also the reference density po 
that will be used for all the calculations in the remainder 
of this paper. Compared with the previously mentioned 
result from molecular dynamics simulations, this result 
is an underestimation of the freezing density by approx- 
imately 7 %. This difference could be due to multiple 
sources of error, including approximating F xs by Eq. ([H]), 
the PY closure relation or not including the perfect lat- 
tice constraint. However, for our purposes, the result is 
acceptable, especially given that the width of the coex- 
istence gap Ap* = (p s — pi) 1 1 pi ss 2% [5l[ is reproduced 
well: Our DFT result is Ap* = 2.20%. 
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From the direct correlation function at freezing point, 
we find the PFC model parameters k m — 6.3965/ct, 
C(k m ) = 0.7855, E s = 14.5487 and E B = 63.7814. The 
direct correlation function, together with the expansions 
of EOF (Eq. (H3J|) and FOF (Eq. $n§ with E B = 0) are 
shown in Fig. [2j As mentioned in the previous section, 
for the EOF and FOF models, the parameters a and b 
are then defined such that the solid phase coexists with 
the liquid phase at p = po, and the amplitude of Fourier 
modes corresponding to the first star of reciprocal lattice 
vectors of the solid phase 

u s = J drn(r)e iG r , (28) 

where G is any reciprocal lattice vector from the first 
star, equals that obtained from the DFT, u s = 0.7914. 
Based on numerical iteration, these two constraints yield 
a = 0.8082 and b = 0.1388 for the EOF, which are no- 
tably rather close to the values of and bid presented 
in the previous section. For FOF, the same fitting proce- 
dure results in only slightly different numbers, a — 0.7812 
and b = 0.1438. For PFC1, the parameter a is chosen to 
fulfill only first of the constraints for EOF and FOF, i.e. 
that the solid coexists with the liquid at p = po. This 
yields a = 1.1934. 

As expected 0, E3] j the solid state density profiles we 
find from the DFT are much more peaked around the 
lattice sites. From all PFC models, we find the field n re- 
sembles the one-mode approximation rather closely. The 
n fields found from EOF and FOF are very similar to each 
other, and to the field obtained from the DFT through 
filtering the higher order Fourier modes in the density. 
The solution from PFC1 differs from the two other PFC 
models in that the amplitude of density fluctuations is 
u s = 0.2051, which is smaller than in the DFT and the 
other PFC models by about a factor four. The coexis- 
tence gap Ap* is 1.57 % in the EOF, 7.70 % in the FOF 
and 0.68 % in the PFC1. Comparing the coexistence gaps 
of the PFC models with the previously mentioned results 
from molecular dynamics and DFT, the EOF gives the 
closest, although not perfect result. In FOF, too small a 
bulk modulus results in too large a coexistence gap [l7j | . 
In PFC1, the small u s and Ap* indicate that the transi- 
tion from solid to liquid is a weaker first order transition 
than in the two other models. 

We have also studied properties of the solid-liquid in- 
terface in these models in the close-packed [10] direction. 
Due to scaling properties of the potential, the interfacial 
free energy of the system is given by 

k B Tfk B T\ 1/6 ~ , . 

r = -2-f-^-J r, (29) 

where T is the dimensionless interfacial free energy in 
rescaled units. Initializing a system with a slab of solid at 
p s in the middle of a liquid at p\ , we find after numerical 
free energy minimization T = 0.234 in the DFT. From the 
EOF and FOF we find almost exactly the same interfacial 



4 r 




-8 -6 -4 -2 2 4 6 8 10 12 

x (atomic layers) 



FIG. 3: (color online) On the top, we show density profiles 
of the solid-liquid interface, from top to down, in the DFT, 
EOF, FOF and PFC1 models. Darker shades of gray corre- 
spond to larger densities in these images, with a scale such 
that the maximum of n in each case corresponds to black, and 
minimum to white. In the lowest image, we show the field n 
averaged in the direction parallel to the interface, as a func- 
tion of the spatial coordinate perpendicular to the interface. 
Black line corresponds to DFT, red dashed line is EOF and 
green line is PFC1 (FOF is not shown, because it would be 
practically indistinguishable from EOF in this plot). 



energies, i.e. T = 0.222 in the EOF and T = 0.223 in 
the FOF. These results also agree rather well with the 
DFT result. On the other hand, from PFC1 we find 
r = 0.0086, which is more than an order of magnitude 
smaller than in all the other models. Density profiles of 
the interface layer in the different models are shown in 
Fig. El It can be seen that in PFC1, the interface layer 
is considerably wider than in all the other models, while 
the interface widths in EOF and FOF are very similar 
to that in the DFT. We are not aware of any computer 
simulation predictions for the surface free energy of r~ 12 
disks. 

As an aside, we note that the justification for a and b 
being different from unity presented in the previous sec- 
tion is not the only one published. Berry et al. [27J have 
noted that a local n 3 term can also be justified by con- 
sidering the k = k' = contribution from a third order 
term in the density expansion of F xs . Technically, this 
is equivalent to assuming that the density field is slowly 
varying compared with the range of three body correla- 
tions. If such term is included in Eq. (j8|), and the loga- 
rithm is expanded in a Taylor series, it is straightforward 
to derive an explicit expression for a, 

a= 1 +p2c( 3 '(0,0), (30) 
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where c( 3 )(0,0) is the k = k' — mode of the three 
body direct correlation function. Using the Ornstein- 
Zernike relation with Percus-Yevick closure, 6^(0, 0) can 
be calculated by noting that it is related to the k = 
mode of the two body direct correlation function through 
the sum rule 



a(3) 



(0,0) 



d6 (2) (Q;po) 



(31) 



The prediction for a we obtain from this approach is 
a w —410, which is not only large in terms of absolute 
value, but also has the wrong sign with respect to stabiliz- 
ing the solid phase. Similar consideration for the fourth- 
order term b, considering k — k' = k" = contribution 
from the four-body correlation term, yields b « 1.35 x 10 3 . 
Based on these considerations, we conclude that our val- 
ues for a and b cannot be justified in terms of higher 
order correlations. 



IV. DYNAMICAL SIMULATIONS 

The solidification front dynamics in the DDFT and the 
different PFC models were studied by growing a hexago- 
nal crystal from an undercooled liquid (i.e. a liquid with 
an initial density pi > pi) in the [10] direction. In the 
direction perpendicular to the solidification front propa- 
gation, the size of the array in our computations is ex- 
actly one interparticle spacing, and periodic boundary 
conditions are used. Due to the periodic boundaries, our 
simulations represent an infinitely wide crystal seed that 
propagates into the liquid. That the size of our simu- 
lation box is only one interparticle spacing in the [01] 
direction means that no instabilities that could roughen 
the surface are allowed. Initial condition for the DDFT, 
EOF and FOF simulations is such that eight monolayers 
of perfect solid are placed in the middle of the under- 
cooled liquid, with a slight smoothing in the boundary of 
solid and liquid phases, while for PFC1, we had to use 
a crystal seed of 12 monolayers in order to initialize the 
growth process at even the smallest undercoolings stud- 
ied. Once the simulation starts, the solid seed grows in 
both directions, and we measure its position as a func- 
tion of time. Position of the surface is defined as the point 
where a line drawn through the local maxima of the den- 
sity corresponding to the solid particles reaches one half 
of its maximum value. In the direction of growth, the 
length of our array was usually 512 interparticle spac- 
ings. 

For numerically integrating the EOF and FOF models, 
we use the well-known semi-implicit operator splitting 
method [28| with Fast Fourier transforms. For spatial 
resolution, we use Ax — \Z3d nn /16 and Ay = d nn /16 
(d nn is the nearest- neighbor distance), time step is At = 
10~ 3 , and the Laplace operator is discretized in fc-space 
as V 2 . = — k 2 . For the DDFT, we employ a similar proce- 
dure, by treating the V 2 n term in Eq. §§§ implicitly and 
the term related to F xs explicitly. The non-local term 



is evaluated in fc-space and derivatives in x and y direc- 
tions are discretized in fc-space as ik x ^ y . For DDFT, in 
order to resolve the sharp density peaks, the linear spatial 
resolution of the PFC is doubled, i.e. Ax = ^/3d nn /32 
and Ay = d ?m /32. The time step we use for DDFT 
is At = 10~ 3 in the regime of low undercooling. In 
the regime of high undercooling we found that retain- 
ing the numerical accuracy required us to decrease the 
time step to At = 10~ 4 , which is smaller than we uti- 
lized in EOF and FOF models by an order of magni- 
tude. This, together with the difference in spatial reso- 
lution, means that using our methods, simulations with 
the DDFT are approximately two orders of magnitude 
slower than with the EOF and FOF models. For PFC1, 
we modified the method used for DDFT such that the 
implicitly treated part is V 2 (l — aCpFci(^ 2 )) n , leaving 
aV ■ nW(CpF ci(V 2 )n) treated explicitly. Even though 
for PFC1 this modification brought great advantage in 
numerical stability, handling the non-linear part still in- 
volves explicit evaluation of sixth derivative. Therefore, 
we found the PFC1 to be numerically most unstable 
among the models studied. In order to ease the require- 
ment this model places on the time step, we dropped 
the spatial resolution perpendicular to growth below that 
used in other PFC models, to Ay = d nn /8, which we did 
not find to have any profound effect on any results pre- 
dicted by the model. However, this only allowed us to uti- 
lize time steps that are one order of magnitude smaller 
than in the DDFT, making the PFC1 numerically ap- 
proximately equally demanding to DDFT using current 
methods. In addition to the differences in Aa;, Ay and 
At between the models, we note that progressing a sin- 
gle time step in DDFT and PFC1 models requires a total 
of five Fourier (or inverse Fourier) transforms, where in 
EOF and FOF, only three are required. 

If the undercooling is small, such that pi < p s , then 
formation velocity of the solid (whose density is always 
at least p s to be stable) is expected to be limited by 
the diffusion of mass to the interface, in an analogous 
manner to the more commonly considered case where 
growth of the solid is limited by transport of heat away 
from the surface . Density of the solid seed in these 
simulations is chosen to be that of the solid coexisting 
with the liquid. As the planar solidification front propa- 
gates, the layer through which diffusion must take place 
widens, and therefore one expects the solidification front 
to propagate as x ~ t 5 . Diffusion controlled growth in 
the PFC model has been previously studied by Tegze et 
al. [6|, but to our knowledge, no studies of the subject 
utilizing a non-local DDFT have been published. On the 
other hand, when density of the liquid from which the 
solid is formed exceeds p s , propagation of the solidifica- 
tion front does not require diffusion of additional mass 
to the surface. Therefore, one expects the front to prop- 
agate with a constant velocity, i.e. x ~ t, that depends 
on the attachment rate of particles on the interface. 

The different regimes of growth, as well as differences 
between the different models, are illustrated in Figs. 0] 
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FIG. 4: (color online) Interface position as a function of time 
when pi — 1.0044/9;. Black solid line is the result obtained 
from the DDFT, red dashed line is from the EOF, blue dash- 
dotted line is from the FOF and green solid line is from the 
PFC1. 



FIG. 6: Growth exponents obtained from the different mod- 
els. Circles are results from DDFT model, squares from EOF, 
diamonds from FOF and triangles from PFC1, while dashed 
line shows the ideal behavior. 
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FIG. 5: Interface position as a function of time when pi 
1.088p;. Different lines are as in Fig. [4] 



FIG. 7: (color online) Logarithmic plot of rescaled interface 
positions as a function of rescaled time from different models 
in A < 0.6 regime. The black dashed line shows the expected 
slope, while other lines are labeled as in Fig. [4] 



and [5] In Fig. [5] we show the interface position as 
a function of time from all models, for a case where 
Pi = 1.0044/?/, which is in the A < 1 regime of all mod- 
els. At the very beginning of the solidification process, 
kinetics define how fast the particles attach to the inter- 
face from the liquid with a density pi. During the initial 
stages, a depletion layer is formed in front of the liquid. 
As time goes on, width of the depletion layer increases, 
and x{t) approaches the expected x ~ t 5 behavior. By 
close inspection of Fig. [4] it can be seen that at the very 
beginning of the process, interface motion in the PFC1 
model is slightly faster than in the DDFT. On the other 
hand in the EOF and FOF models, the initial interface 
motion is slightly slower than in DDFT, indicating more 
restriction to growth due to interface kinetics. In the 



EOF and FOF models, the initial velocities are strik- 
ingly similar. On the other hand, Fig. [5] shows interface 
positions as a function of time for an initial density of 
Pi = 1.088/0;, which is in the A > 1 regime of all models. 
Again, in the initial stages of solidification, a depletion 
layer is formed in front of the moving interface. However 
in this case the width of the depletion layer, and thus the 
propagation velocity of the interface, quickly approach 
constant values, and therefore the growth seems linear. 

In order to quantify the detection of the different 
growth regimes, we have fitted a power law growth func- 
tion, 

x = x + c(T-T ) a , (32) 
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FIG. 8: Growth rates d obtained in the diffusion controlled 
regime. Circles are results from DDFT model, squares from 
EOF, diamonds from FOF and triangles from PFC1, with 
lines connecting the symbols. Inset shows scaled data. 



in the surface positions as a function time resulting from 
the different models at different initial densities. The first 
quarter of the x(t) data is ignored in these fits in order 
to minimize the effect from the initial stages, while still 
obtaining a fairly robust fit for the four free parameters. 
The exponents a resulting from these fits are shown in 
Fig. HI The growth exponents obtained from the FOF 
and PFC1 models are in general close to the ideal re- 
sults (even at unit undercooling, the observed exponent 
a Ri 0.70 from FOF is close to the expected anomalous 
exponent a = 2/3 [H, Hl)> whereas in DDFT and EOF 
models, the transition from the diffusion controlled to 
kinetics controlled regime seems more continuous. The 
discrepancy between observed and expected exponents is 
due to insufficient simulation time for the formation of 
the steady-state depletion layer in DDFT and EOF mod- 
els. The discrepancy is most evident in the EOF model. 
This is most likely because the combination of small co- 
existence gap, fast diffusion in the solid, and relatively 
slow interfacial kinetics makes the formation of a quasi- 
steady-state depletion layer in the EOF slowest among 
the models. 

In the A < 1 regime, we have quantified the effect of 
the initial density on the quasi-steady-state velocity of 
front propagation by fitting the x(t)'s resulting from dif- 
ferent models for different p^s with the expected growth 
law, 



x + d(r - t ); 



(33) 



ignoring the early stages where xa~l < 50. Eq. (f3"3"]) is 
mathematically equivalent to Eq. (|3"2"|) with a set to 0.5. 
Despite the previously mentioned discrepancies between 
the observed and now pre-set growth exponents, we are 
able to obtain good fit with Eq. (|3"3")l for all models in 
the regime A < 0.6, as shown in Fig. [7J The d resulting 



FIG. 9: Growth velocities v obtained in the kinetics controlled 
regime. Symbols as in Fig. [8] 



from these fits are shown in Fig. [5] Beyond A = 0.6, we 
have less confidence in having reached close enough to 
a steady-state diffusion-controlled growth regime, espe- 
cially in the EOF model, and therefore the data is only 
shown for undercoolings up to A = 0.6. 

It is seen in Fig. [8] that among the PFC models, 
the EOF in general gives the closest agreement with the 
DDFT. Due to scaling properties of the problem, if the 
growth were purely diffusion limited, one would expect 
the position of the interface to depend on the dimen- 
sionless undercooling A and \/ Dt, where D = 1 — (5(0) 
(for PFC1, multiply 6(0) by a) is the effective diffu- 
sion constant in a given model (in the limit of small, 
long- wavelength density fluctuations, with this definition, 
all the models studied reduce to the diffusion equation 
d T n = DV 2 n). Such a scaling law for the growth rates is 
illustrated in the inset of Fig. [5J where we show that the 
cf s scaled by VD as a function of A in DDFT, EOF and 
FOF models follow the same curve, indicating that dif- 
ferences in microscopic details of those models are unim- 
portant in determining the front velocity in the diffusion 
controlled regime. Results from PFC1 lie slightly above 
those from the other models in the rescaled plot, which 
we believe is most likely a result of numerical inaccuracy. 
These scaling properties of the problem are the reason 
why the EOF, which has exactly the same D and close 
to the same Ap* as DDFT, reproduces the result of the 
DDFT with higher accuracy than the FOF, which re- 
sults in a d that is approximately an eighth of the result 
obtained from the DDFT for the same pi. The scaling 
argument also suggests that the close agreement of PFC 1 
to EOF and DDFT in the small density limit probably 
results from a cancelation of errors due to the smaller D 
and smaller Ap* in the PFC1 model. 

In the regime where A > 1 , the results have been fitted 
with a linear growth law, 



x — Xq + VT, 



(34) 
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where the effect of the initial transient was removed by 
ignoring data for which x < 100c. The resulting v as a 
function of pi from all the models are shown in Fig. GO As 
expected, the front velocity increases as the initial density 
is increased in all the models. It is also apparent that for 
any given initial density, the velocity obtained from the 
EOF is the closest approximation to the DDFT among 
the present PFC models. If the density axis is rescaled by 
subtracting the density of the solid coexisting with the 
liquid, the results from FOF seem to agree with DDFT 
practically as well as those from the EOF, as shown in the 
inset of Fig. [§] On the other hand, the velocities observed 
in PFC1 model seem to be a significant overestimation 
when compared with the results from all the other models 
studied, even after rescaling the densities. This suggests 
that the smaller amplitude of density fluctuations and 
wider crystal-melt surface result in a kinetic barrier which 
is somewhat smaller than in the other models. 



V. CONCLUSIONS 

We have presented a new way to derive the eighth- 
order phase field crystal model (EOF) from the density 
functional theory of classical systems. The model was 
applied to study solidification front dynamics in a two- 
dimensional ensemble of particles interacting via r~ 12 
potential. Predictions from the EOF were compared 
with similar predictions from dynamical density func- 
tional theory (DDFT) of Marconi and Tarazona, and two 
previously presented phase-field crystal (PFC) models. 
For the static properties of the system in these mod- 
els, we find that the DFT predicts freezing of the r~ 12 
disks at a density that is about 7 % lower than seen in 
molecular dynamics simulations. From the PFC models 
studied, we find that the EOF gives the most accurate 



description of the static properties of the material under 
study. By studying crystal growth in the diffusion con- 
trolled regime, we find that the EOF gives the best agree- 
ment with DDFT among the phase field crystal models, 
due to the most accurate description of liquid diffusion 
constant and solid-liquid coexistence gap in the model. 
In the regime of interface kinetics controlled growth, we 
again find the EOF gives closest agreement to the DDFT 
for all initial densities, although if the initial density is 
rescaled by the melting point of the solid, the fourth or- 
der fitting scheme slightly outperforms the EOF. These 
results suggest that among the PFC models studied, the 
EOF gives the closest approximation to the DDFT. This 
implies that the EOF is a good candidate for a model 
to be used for atomistic scale simulations of the growth 
of two dimensional hexagonal crystals of Brownian par- 
ticles, at least in the absence of an external field. In the 
presence of an external field, a further study would be 
required to quantify the response in the different models. 
It should also be noted that while the current study has 
considered a simple two-dimensional problem, a similar 
study of the growth of a three dimensional crystal could 
also quantify the differences in anisotropy of the different 
phenomena in the models. 
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